clear all
set more off

// Robustness test: Capitals close to dep borders

global input "...\replication\data"
global output "...\replication\outputs"

use "$input\com_data.dta", clear


****************** Costmetics ******************
** Reomove Corsica
ren com_res com
drop if substr(com,2,1) == "A"
drop if substr(com,2,1) == "B"
destring com, replace

****************** Get dist border ******************
merge m:1 com using "$input\com_arcgis.dta"
drop if _merge != 3
drop _merge
** Drop those within 30km of coast or national border
keep if dist_nat_border < 0 


****************** Get dist instrumental border ******************
merge m:1 com using "$input\com_dist_instr.dta"
drop if _merge != 3
drop _merge
merge m:1 com using "$input\same_dir.dta"
drop if _merge != 3
drop _merge

***** Focus on those with close to border capital ****
** 17 (Charente Maritime), 79 (Deux Sevres), 65 (Hautes Pyrenees), 84 (Vaucluse), 03 (Allier), 71 (Saone et Loire), 58 (Nievre), 57 (Moselle), 26 (Drome)

gen rob_gp = 0
foreach i in 17 79 65 84 03 71 57 58 26 {
replace rob_gp = 1 if dep_res == "`i'"
replace rob_gp = 1 if dep_work == "`i'"
}
keep if rob_gp == 1

****************** Get main vars******************
** Outcome
egen sum_ind = sum(ind), by(com)
gen p_ind = ind/sum_ind
** Treatment
gen Treat = 0
replace Treat = 1 if same == -1

drop if dist_border == 0

keep if abs(dist_border) <= 25

gen dist_instr2 = dist_instr
replace dist_instr2 = -dist_instr2 if dist_border < 0
replace dist_instr2 = -dist_instr2 if same_dir == 0

//////////////// Count ///////////////////////
destring dep_res dep_work, replace


****************** calonico et al ******************
rdbwselect p_ind dist_border, c(0) 
local opt_h=e(h_mserd)

rdrobust p_ind dist_border ,c(0) h(`opt_h') 
sca N_bw=e(N_h_l)+e(N_h_r)
estadd scalar N_bw=N_bw
sum p_ind if abs(dist_border) < `opt_h' & dist_border < 0
sca m_outc = r(mean)
estadd scalar m_outc=m_outc
est store ATE_1

rdrobust p_ind dist_border ,c(0) h(`opt_h') covs(dep_res dep_work) 
sca N_bw=e(N_h_l)+e(N_h_r)
estadd scalar N_bw=N_bw
sum p_ind if abs(dist_border) < `opt_h' & dist_border < 0
sca m_outc = r(mean)
estadd scalar m_outc=m_outc
est store ATE_2

rdrobust p_ind dist_border ,c(0) h(2*`opt_h') 
sca N_bw=e(N_h_l)+e(N_h_r)
estadd scalar N_bw=N_bw
sum p_ind if abs(dist_border) < 2*`opt_h' & dist_border < 0
sca m_outc = r(mean)
estadd scalar m_outc=m_outc
est store ATE_3

rdrobust p_ind dist_border ,c(0) h(0.5*`opt_h') 
sca N_bw=e(N_h_l)+e(N_h_r)
estadd scalar N_bw=N_bw
sum p_ind if abs(dist_border) < 0.5*`opt_h' & dist_border < 0
sca m_outc = r(mean)
estadd scalar m_outc=m_outc
est store ATE_4

****************** Results in Table ******************
esttab ATE_1 ATE_2 ATE_3 ATE_4, tex ///
b(%9.3f) se(%9.3f) stats( h_l N N_bw m_outc, labels("Bandwidth" "Obs" "Obs. bw." "Mean dep. var. below" ) ///
fmt( %9.3gc %9.0gc %9.0gc %9.3gc)) nolines nogaps star(* 0.10 ** 0.05 *** 0.01) ///
keep(RD_Estimate) coef(RD\_Estimate  "ATE") 


***********************************
************************** IV ***********************************
***********************************
rdrobust p_ind dist_border ,c(0) h(`opt_h') fuzzy(dist_instr2)
sca N_bw=e(N_h_l)+e(N_h_r)
estadd scalar N_bw=N_bw
sum p_ind if abs(dist_border) < `opt_h' & dist_border < 0
sca m_outc = r(mean)
estadd scalar m_outc=m_outc
est store ATE_1

rdrobust p_ind dist_border ,c(0) h(`opt_h') covs(dep_res dep_work) fuzzy(dist_instr2)
sca N_bw=e(N_h_l)+e(N_h_r)
estadd scalar N_bw=N_bw
sum p_ind if abs(dist_border) < `opt_h' & dist_border < 0
sca m_outc = r(mean)
estadd scalar m_outc=m_outc
est store ATE_2

rdrobust p_ind dist_border ,c(0) h(2*`opt_h') fuzzy(dist_instr2)
sca N_bw=e(N_h_l)+e(N_h_r)
estadd scalar N_bw=N_bw
sum p_ind if abs(dist_border) < 2*`opt_h' & dist_border < 0
sca m_outc = r(mean)
estadd scalar m_outc=m_outc
est store ATE_3

rdrobust p_ind dist_border ,c(0) h(0.5*`opt_h') fuzzy(dist_instr2)
sca N_bw=e(N_h_l)+e(N_h_r)
estadd scalar N_bw=N_bw
sum p_ind if abs(dist_border) < 0.5*`opt_h' & dist_border < 0
sca m_outc = r(mean)
estadd scalar m_outc=m_outc
est store ATE_4

****************** Results in Table ******************
esttab ATE_1 ATE_2 ATE_3 ATE_4, tex ///
b(%9.3f) se(%9.3f) stats( h_l N N_bw m_outc, labels("Bandwidth" "Obs" "Obs. bw." "Mean dep. var. below" ) ///
fmt( %9.3gc %9.0gc %9.0gc %9.3gc)) nolines nogaps star(* 0.10 ** 0.05 *** 0.01) ///
keep(RD_Estimate) coef(RD\_Estimate  "ATE") 

